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Abstract. The joint spectral radius of a set of matrices is a measure of the maximal 
asymptotic growth rate that can be obtained by forming long products of matrices taken 
from the set. This quantity appears in a number of application contexts but is notoriously 
difficult to compute and to approximate. We introduce in this paper a procedure for 
approximating the joint spectral radius of a finite set of matrices with arbitrary high 
accuracy. Our approximation procedure is polynomial in the size of the matrices once the 
number of matrices and the desired accuracy are fixed. 

For the special case of matrices with non-negative entries we give elementary proofs of 
simple inequalities that we then use to obtain approximations of arbitrary high accuracy. 
From these inequalities it follows that the spectral radius of matrices with non-negative 
entries is given by the simple expression 

p(Ai, . . . , A^) = hm p'l\Af + • • ■ + ) 

fe— >oo 

where it is somewhat surprising to notice that the right hand side does not directly involve 
any mixed product between the matrices {^A®^ denotes the fc-th Kronecker power of A). 

For matrices with arbitrary entries (not necessarily non-negative) we introduce an ap- 
proximation procedure based on semi-definite liftings that can be implemented in a recur- 
sive way. For two matrices, even the first step of the procedure gives an approximation 
whose relative accuracy is at least that is, more than 70%. The subsequent steps 

improve the accuracy but also increase the dimension of the auxiliary problems from which 
the approximation can be found. 

Our approximation procedures provide approximations of relative accuracy 1 — e in time 
polynomial in j7,(iii ^ where m is the number of matrices and n is their size. These bounds 
are close from optimality since we show that, unless P=NP, no approximation algorithm 
is possible that provides a relative accuracy of 1 — e and runs in time polynomial in n and 
1/e. 

As a by-product of our results we prove that a widely used approximation of the joint 
spectral radius based on common quadratic Lyapunov functions (or on ellipsoid norms) 
has relative accuracy where m is the number of matrices.^ 
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initiated by the Belgian State, Prime Minister's Ofhce, Science Policy Programming. The scientific 
responsibility is assumed by the authors. 
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1. Introduction 



Let II • II be a matrix norm. The spectral radius of the real matrix A is defined by 



The spectral radius of a matrix does not depend on the chosen matrix norm and it is 
also equal to p{A) = max{|A| : A is an eigenvalue of A} (see for example Corollary 
5.6.14]). The definition of spectral radius can be extended to sets of matrices in a natural 
way. The joint spectral radius of a set of matrices is a quantity introduced by Rota and 
Strang in the early 60 's that measures the maximal asymptotic growth rate that can be 
obtained by forming long products of matrices; see The general definition is for 

arbitrary sets of matrices but we shall consider here only finite sets. Let {Ai, . . . , Am} be 
some set of real matrices. To the finite sequence a = (cxi, (T2, . . . , cxfc) e {1, . . . , m}^ we 
associate the corresponding matrix product 

A — A ... 4 A 

With this notation, the joint spectral radius is defined by 



As for the single-matrix case, the joint spectral radius does not depend on the matrix 
norm used. To see this, remember that any two matrix norms ||.||i and ||.||2 are related 
by a||A||i < ll^lh < /3||^||i for some < a < (3. For any product a G {1, . . . ,m}^ one 

has a-'^/'^IIAo-lli'''' < ||^o-||2^'^ ^ /^"^'''^ II ^cr II i^'^ and by letting k tend to infinity we conclude 
that the joint spectral radius is well defined independently of the matrix norm used. 

A definition analogous to ()1.2j) is possible by replacing the norm appearing in the def- 
inition by a spectral radius. The quantity defined in this way is the generalized spectral 
radius introduced in In |3j, the joint and generalized spectral radii of finite (or 

bounded) sets of matrices are proved to be equal (see also |16j Theorem 1] for an ele- 
mentary proof); thus reinforcing the status of the joint spectral radius as a legitimate 
generalization of the spectral radius of a single matrix. In the sequel we shall deal only 
with the spectral radius defined with a norm, as in ()1.2|) . 

Since its introduction in the 60's the joint spectral radius has appeared in a number 
of different contexts; see, e.g., [HIl or [12] for recent short surveys^. Let us illustrate 
one application in a dynamical system context. Consider the simple discrete-time linear 
inclusion 

Xk+lE {Ai,...,Ara} Xk, Xo G R" 

in which at each step a particular linear transformation is chosen from of a finite number 
of possible choices. The maximal asymptotic rate of growth of the trajectories associated 
to such a system is given by the joint spectral radius p{Ai, . . . ,Am)- (In the discrete 
linear inclusion literature, the logarithm of the joint spectral radius is sometimes called 

^Google returns 625 entries upon entry of the query "joint spectral radius" . 



(LI) 



p{A) = lim WA'^W^/''. 



(1-2) 



p{Ai, . . . , Am) = lim sup max 
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Lyapunov indicator, see for example |2j.) In particular, all possible trajectories will con- 
verge to the origin if and only if p(Ai, . . . , Am) < 1- The joint spectral radius can thus be 
associated with the stability properties of time-varying linear systems in the worst case 
over all possible time variations. It also occur in the context of "asynchronous" jSS] or 
"desynchronised" [211 systems. Besides systems-theoretic interpretations, the concept is 
pervasive in many areas of applied mathematics such as in wavelets ^3] , iterated function 
systems, random walks, fractals, numerical solutions to ordinary differential equations 
[T5] . discrete-event systems interpolation [HH], and coding theory j2Hl- 

Despite its natural interpretation, the joint spectral radius is difficult to compute. 
Questions related to its computability and to the existence of efficient approximation 
algorithms have been posed more than a decade ago (see ES])- In principle, the 

spectral radius can be approximated to any desired accuracy by computing converging 
sequences of upper and lower bounds. The following bounds, proved in (25j, 



can be evaluated for increasing values of k and lead to arbitrary accurate approximations 
of p (see, e.g., [1^] or [I7j). Such approximation algorithms can in turn be used in proce- 
dures that decide, after finitely many steps, whether p > 1 or p < 1 (such procedures are 
given, e.g., by Brayton and Tong JH] in a system theory context and by Barabanov ^ in 
the context of discrete linear inclusions). These procedures may however not terminate 
when p happens to be equal to 1 and the existence of algorithms for computing arbitrarily 
precise approximations of p does therefore not rule out the possibility that the decision 
problem "p < 1?" is undecidable. It is so far unknown whether this is the case or not (see 
||25j for a discussion of this issue and for a description of its connection with the finiteness 
conjecture that has since then been proved to be false; see \^A^ and [4,). A negative result 
in this direction is given in jH] where it is proved that the related problem "p < 1?" is 
algorithmically undecidable. 

Approximations of the joint spectral radius that are directly based on the inequalities 
fll.3|) are expensive to compute. In [^Zl, the exponential number of products that appear 
in the direct and naive computation of the bounds in p.3|) is reduced by avoiding dupli- 
cate computation of cyclic permutations; the total number of product to consider remains 
however exponential. In addition to this, there are to this date no known theoretical 
guarantees for the rate of convergence of the bounds appearing in (jl.3|) . Approximations 
of arbitrary degree of accuracy can be computed, but at a price that may happen to be 
prohibitive. 

There are in fact intrinsic limitations for the rate at which the joint spectral radius 
can be approximated. Let us say that the value ^ approximates the value p with relative 
accuracy p G [0, 1] (or 100 p %), if p ^ < p < ^. By using a small adaptation of a proof 
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appearing in |32] we show that, unless P=NP, there is no algorithm that can compute 
the joint spectral radius of two matrices with relative accuracy 1 — e in time polynomial 
in the size of the matrices and in 1/e (see later for more precise definitions). Despite this 
negative result it is still conceivable, as pointed in [Hj, that for any fixed desired accuracy, 
there exists a polynomial time algorithm that computes the joint spectral radius of the 
matrices with that accuracy. We prove in this paper that this is indeed the case. In Theo- 
remEl we prove that the joint spectral radius of m matrices of size n can be approximated 
with relative accuracy 1 / i/m by computing the spectral radius of a single matrix whose 
size is less than n^. This procedure can be applied in a recursive way and in general we 
show how a relative accuracy of {l/^/my/'' can be obtained by computing the spectral 
radius of a single matrix of size less than n^'^. As an illustration, the spectral radius of 
two matrices of size n can be computed with an accuracy of 70% by computing the spec- 
tral radius of a single matrix of size n^, and an accuracy of 95% can be obtained for the 
joint spectral radius of three matrices by computing the spectral radius of a single matrix 
of size n^^. More generally, for any number of matrices in the set and desired relative 
accuracy we construct in Section |3] a single matrix whose spectral radius approximates 
the joint spectral radius with the desired accuracy. The approximation procedure runs 
in polynomial time once the desired accuracy and the number of matrices in the set are 
fixed. More precisely, our approximation procedures provide approximations of relative 
accuracy 1 — e in time polynomial in n*^''^™-'/^, where m is the number of matrices and n 
is their size. Notice that n*^'"^"^^/^ = e(inninm)/€ _ ^{inn)/e ^^^^ approximation pro- 

cedure also runs in time polynomial in the number of matrices once the desired accuracy 
and the size of the matrices are fixed. These bounds are close from optimality since we 
prove that, unless P=NP, no approximation polynomial time algorithm in n and 1/e is 
possible. 

We now briefly describe how our results are obtained and how the paper is organized. 
We first notice in Section |21 that for matrices with non- negative entries the joint spectral 
radius satisfies 

(1.4) ^ p{Ai + --- + Am)< p{Ai, ...,Am)< p{Ai + --- + Am). 

(The left-hand side inequality is valid for arbitrary matrices; there is no need to assume 
that the matrices have non- negative entries.) Matrices with non- negative entries are 
exactly those matrices A that are such that A R " C R " and we prove in Theorem ^ 
that the inequalities (jl.4|) are satisfied not only for sets of matrices with non-negative 
entries but also for sets of matrices that leave a proper cone invariant, i.e., matrices Ai 
that are such that Ai K C K for some proper cone K and all i. Thus for these matrices 
the inequalities p.4|) provide a relative accuracy 1/m. This accuracy can be improved by 
considering Kronecker powers of matrices. In Section El we give an elementary proof that 
the joint spectral radius of the k-th Kronecker powers of the matrices in a set is equal to 
the k-th power of the joint spectral radius of the set (Theorem E)). Combining this with 
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the relation ()1.4j) we prove that the approximation 



P 



1 



has relative accuracy Some of the arguments used to derive this result are valid 

only for matrices that have a common proper invariant cone. In Section 0] we show how 
similar arguments can be used for arbitrary matrices. We introduce a semi-definite lifting 
procedure that transforms a linear operator acting R" into a linear operator acting on 
the space of symmetric n x n matrices. It appears that under this transformation the 
joint spectral radius is simply squared. Moreover, the lifting has the interesting feature 
that the operator defined in this way leaves the cone of semi-definite matrices invariant. 
This observation leads to a simple approximation procedure (Theorem El) for arbitrary 
matrices. Implementation issues and numerical examples are provided in Sectional Our 
approximation procedure provides a relative accuracy 1 — e in time polynomial in jt.C'i™)/'^. 
We prove in Section El (Theorem IH} that the same accuracy cannot be obtained in time 
polynomial in n and 1/e unless P=NP. The semi-definite lifting introduced in Section HI 
is a fundamental and powerful tool for stability analysis. As an illustration of this, we 
prove in Section [71 how to apply our results to analyze the quality of the so-called elhpsoid 
approximation of the joint spectral radius. We prove that the ellipsoid approximation (a 
notion that is formally defined in jlj but that is implicitly present in a number of earlier 
contributions on hybrid and time- varying systems, see, e.g., j2I] and as well as some of 
the references in [26j) is an approximation of guaranteed accuracy 1/y/m; a proof of this 
corollary can also be extracted from Section 3 of PTj. For matrices of high dimension this 
result significantly improves the earlier bound of 1/ ^/n proved in [4] and in ^ . Finally, 
in a last section, we discuss some of our results. 



In this section, we consider sets of matrices that leave a proper cone invariant and we 
show with elementary arguments how joint spectral radius approximations of guaranteed 
accuracy can easily be computed for this case. We start with a proof that the spectral 
radius of a convex combination of matrices is always less or equal to the joint spectral 
radius of the matrices. This result is valid for all sets of matrices (there is no need to 
assume that the matrices have non-negative entries) and is proved in JH] using the main 
result of [241. Here we present a direct and elementary justification. 

Lemma 1. For any set of matrices {Ai : i = l,...,m} and any > satisfying 
Yli=i ai = 1 we have: 



2. Approximation for matrices leaving a cone invariant 



m 



(2.1) 
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Proof. Let us fix some > with Yli Oii = 1 and an integer k>l. Then 

m 

\\{S^aiAif\\ = \\ V a„A„\\< V a,\\A,\\ < max \\A,\\ 

•^"^ •^"^ (Te|l,...,m}* 

«=1 (TG{l,.--,m}'° o-ell,...,™}* 

but then also 

m m 

hm IK^^ajAj)''!!^/*^ = hmsup IKy^ajAj)*^!!^/'' < hmsup max 11^4^1^^'' 
1=1 1=1 

and the result then follows from the definitions of the spectral radius (jl.lll and of the 
joint spectral radius p.2j) . 

□ 

An immediate corollary is given by: 
Corollary 1. 

m 

^ p(^A,) <p(Ai,...,/lJ. 
i=l 

□ 

The example Ai = A,A2 = —A clearly shows that we cannot hope in general to have 

m 

p{Ai, . . . , Ajn) < p(X] ^i)- This inequality is nevertheless satisfied when the matrices A^ 

i=l 

leave a proper cone invariant. A cone in R" is a subset C R" such that \v E K for all 
A > and v (z K. We say that a cone K is proper if it is closed, convex, has nonempty 
interior, and contains no straight line. For example, the set of vectors with non-negative 
entries R" is a proper cone but R" itself is not (see for more background on cones 
and proper cones). We shall say that the matrices A^ leave a proper cone invariant if 
there exists a proper cone ii" C R" such that AiK <ZK for all i. For example, matrices 
with non-negative entries leave the proper cone R" invariant. 

Lemma 2. Associated to any proper cone K there is a matrix norm || ■ that satisfies 
^ + for all matrices A and B that leave the cone K invariant. 

Proof. Most usual matrix norms satisfy this property when K = R"; we provide a 
construction for arbitrary proper cones K. 

Let II ■ II be some arbitrary vector norm. The dual vector norm || ■ ||^, is defined by 
||w||* = max||j,||=i w'^v. Assume that 7^ is a proper cone such that Ai K O K. Define the 
dual of K by 

i^* = {w e R" : w'^v > for all v G K}. 

The dual of a proper cone is again a proper cone, see Let us now consider the 

quantity 

II A II X = max nF Av. 

vGK,m£K* 
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It is easy to verify that because K and its dual are proper cones, the quantity || ■ is 
indeed a matrix norm. Moreover, by the definition of the norm, matrices A and B that 
satisfy AK C K and BK C K, are such that \\A\\k <\\A + B\\k. □ 
We may now state the result of this section. 

Theorem 1. Let {Ai : i = l,...,m} be a set of matrices that leave a proper cone 
invariant. Then 

m m 

(2.2) ^p(E A) < p{Au ...,AJ< piZ Ai). 

1=1 i=l 

Proof. The lower bound in (|2.2j) is already established in Corollary ^ Assume that the 
matrices in {Ai : i = 1, . . . ,m} leave the proper cone K invariant and define the matrix 
norm \\-\\k as in Lemma|21 Since AiK C K, we also have A„K C K for all cr G {1, . . . , m}*^ 
but then also 

max \\AJk < II J2 ^-Wk = 11(52 

(Te|l,...,m|''" ■^■"^ ■^"^ 

a-e{l,...,m}'' 4=1 

and therefore 

m m 

p(Ai, . . . , Am) = limsup max ||A^||]/^ < limsup || (^^ A»)^||)^^ = pC^^Aj). 

k^oo o-g{1, .••,»"}''' fc->oo . -, . 1 

1=1 1=1 

□ 

It is interesting to notice that the conclusion of Theorem ^ does not directly involve 
the cone K. The existence of a proper cone satisfying the hypotheses suffices to conclude 
and the exact nature of the cone is irrelevant. 

A matrix with non- negative entries leaves the cone invariant and so we have the 
following corollary. 

Corollary 2. Let the matrices Ai, i = 1, . . . ,m, have non-negative entries. Then 

m m 

i PiY, A) < • • • > < piY. ^*)- 

1=1 i=l 

□ 



3. KRONECKER LIFTING 

We now describe a way to exploit Theorem^ for obtaining approximations of arbitrary 
high accuracy. Our approximations involve Kronecker powers of matrices. The Kronecker 
product of two matrices A G Rp^^'^^ and B G Rf^^i'^ is the pip2 x qiq2 matrix defined by 

/ ... A,,,,B \ 

A(^B= I 

\ Ap^^iB . . . Ap^^q^B J 
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We also define the Kroneker power of a matrix: 

A^'' = A®A...A®A. 
k times 

Tliere is no need for parentliesis in tliis expression since tfie Kronecker product of two 
matrices is an associative operation. Let us prove some elementary properties. 

Lemma 3. For matrices of appropriate sizes we have: 

(1) {A(^Bf = A^^B^; 

(2) (Ai ® 50(^2 ® B2) = {A1A2) ® (^i^s); 

(3) {A^'^f = 

(4) A'^'^B'^^ = 

(5) Let II ■ II denote the spectral matrix norm induced by the standard Euclidean vector 
norm. Then ||A®*^|| = ||Af\ 

Proof. The first statement directly follows from the definition of the Kronecker product. 
The second statement is quite standard; see, e.g., Section 4.2 in [221 • The third and 
fourth statement follow, respectively, from repeated application of the first and second 
statement. Finally, for proving the last statement, note that by using (3) and (4) we 
obtain 



WA'^'^f = maxx^(A®^)^A®^x = max = max x"" {A^ Af' x . 

||x'j|=l ||x||=l |ja;|j=l 

The matrix (A'^A)'^'' is symmetric and so the last expression is also equal to the largest 
magnitude of the eigenvalues of (A'^A)'^'^. For a matrix B the eigenvalues of the matrix 
B^'' are given by all possible products of k eigenvalues of B (see Theorem 4.2.12 in E2]) 
and we therefore conclude max x'^{A'^A)'^^x = II^H*'. □ 

||x||=l 

We can now prove a useful identity for the spectral radius of Kronecker powers of 
matrices. 

Theorem 2. Let {Ai : i = 1, . . . , m} be a set of matrices and I > 1. Then 

p«,...,Af)=p'(Ai,...,AJ. 
Proof. Let cr G {1, . . . ,m}^ By using the equalities (3) and (5) of Lemma El we obtain 

lK---<ll = ll^fll = 11^.11'- 

The result then follows from the definition of the joint spectral radius (jl.2|) . □ 
We now would like to combine Theorem |21 with Theorem Q in order to obtain approxi- 
mations of arbitrary high accuracy for matrices leaving a cone invariant. In order to do 
this we need to consider Kronecker products of cones. For two sets Qi C i?" and Q2 C R"^ 
let us denote 

Qi<S)Q2 = {z = u®v : u e Qi, V G Q2}- 
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(In this definition, the vectors u and v are treated as column matrices.) If Ki and K2 
are proper cones then Ki K2 does not need to be a proper cone. Indeed, consider for 
example the proper cones Ki = K2 = Rl for which it is easy to see that Ki K2 has 
empty interior in R^. We do however have the following result. 

Lemma 4. Let Ki and K2 he proper cones, then the cone K = Conv {Ki ® K2) is a 
proper cone also. 

Proof. The cone K is closed and convex by definition. Assume that int K = ^. Then 
there exists A G \ {0} such that 

{A,z) =0 yzeKD Ki(^ K2. 

Note that the function f{u,v) = {A,u ® f ) is a bilinear form defined on Ki x K2. Since 
int {Ki X K2) 7^ 0, / can vanish on this set only if A = 0. That is a contradiction. 

Let us prove now that K contains no straight line. Assume that this is not the case. 
Then int K* = 0, which implies existence of 5 G i?"'" \ {0} such that 

(3.3) {B,w)=0 'iweK*. 

Since {u ® v,x ® y) = {u,x) ■ {v,y), we have ^ K2 C K*. Since and have 
nonempty interior, relation (|3.3p is impossible by the first part of the proof. □ 
We are ready to state the main result of this section. 

Theorem 3. Let {Ai : i = 1, . . . , m} be a set of matrices leaving a proper cone invariant. 
Then 

_!_ pi/k^^^k + . . . + < p^A^^ ...,AJ< p'/'{Af + ■■■ + AZ'). 
In particular, the joint spectral radius is given by 

(3.4) p{A^, ...,A^)= hm p'/''{Af + ... + A%''). 

k^oo 

Proof. Assume that the matrices leave the proper cone K invariant. If AK C K, then 
(AK)^'' C K®'' and by Lemma A^'^K®'' C K®''. But then also A'^'' Conv K'^'' C 
Convi^®^. By Lemma 0] the cone Conv K^^ is a proper cone. Thus the matrices Af^ 
leave a proper cone invariant and by Theorem^ we have 

i. p(Af + . . . + ) < p(Af , . . . , ) < p{Af + ■ ■ ■ + ). 

In order to conclude it suffices to use the fact, proved in Theorem that 

p{Af,...,AZ') = p\Ar,...,A^). 

□ 

It is somewhat surprising to notice that the right hand side in ()3.4j) does not directly 
involve any mixed product between the matrices. 

An approximation of relative accuracy (1/m)^/'^ can thus be obtained by computing 
the spectral radius of a single matrix of dimension n^. For pairs of matrices some of the 
relative accuracies and the corresponding matrix sizes are as follows: 
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Relative accuracy 0.707 0.840 0.917 0.957 
Matrix size n^'^ 

4. Semi-definite lifting 

The result in Theorem El provides an easy way to evaluate with arbitrary accuracy 
the joint spectral radius of matrices that leave a proper cone invariant. Unfortunately, 
matrices do not always leave a proper cone invariant. We show in this section how an 
invariant cone can always be created by semi-definite lifting. 

Let A G R"^" and consider the following linear operator: 

(4.1) X ^ AXA^ : R"^" ^ R"^". 

A matrix representation for this linear operator can be obtained by using the matrix-to- 
vector operator that develops a matrix into a vector by taking its columns one by one. This 
operator, denoted vec, satisfies the elementary property \ec{CXD) = {D^ ® C) vec(X) 
(see Lemma 4.3.1 in jlS])- We therefore have Yec{AX A^) = {A® A) vecX and a matrix 
representation of the linear operator (j4.ip is thus simply given hy A® A. 

Consider now the space of symmetric matrices S. This space is a subspace of R"'^" of 
dimension n{n + l)/2 and the operator 

(4.2) X ^ AXA^ : S ^ S. 

is again a linear operator. Matrix representations of this operator are of course different 
from those of the operator ()4.1|) since, in particular, the spaces on which these operators 
operate have different dimensions (we describe in the next section how to construct a 
matrix representations for the operator on symmetric matrices). In the next theorem we 
prove that even though the operators X AXA^ on S and on R,"^"- are different, their 
joint spectral radius are equal. 

Theorem 4. Let {Ai G R"^" : z = 1, . . . , m} and denote by Ma a matrix representation 
of the linear operator X — > AXA^ : S ^ S. Then we have 

piMA,, . . . , MaJ = piAi (g) Ai, . . . , A^ (g) A^) = p\Ai, . . . , AJ. 

Proof. The second equality is already proved in Theorem |21 we only need to prove the 
first equality. Let || ■ \\f denote the Frobenius matrix norm (i.e., the sum of squares of all 
entries) and consider the resulting induced operator norms for X AXA^ on R,"^" and 
S: 

sup II AXA^^IIi? and sup || AXA'^Hi;'. 

xeR"^",||x||_F=i X(^s,\\x\\p=i 

We claim that these two operator norms are equal and that the supremum is achieved 
for some symmetric matrix of rank one. Indeed, using the matrix-to-vector operator vec 



COMPUTATIONALLY EFFICIENT APPROXIMATIONS OF THE JOINT SPECTRAL RADIUS 11 



and denoting by || • || the usual vector Euclidean norm we get 

sup IIAXA'^lll = sup ||(A (g) y4)a;||^ 

= sup x^{AiS) A)'^{Aig) A)x 
= sup x'^iA'^A (g) A'^A)x 

For deriving the last equality, we have used (2) in Lemma 121 The supremum in the 
last expression is achieved by any eigenvector x associated to an eigenvalue of largest 
magnitude of A^A ® A^A. Let (Aj, Vi) i = 1, . . . ,n he the set of eigenvalues/eigenvectors 
pairs associated to the symmetric matrix A^A and let Ai be such that |Ai| > Aj for 
all i. The eigenvalues/eigenvectors pairs of the matrix A^A (g A^A are then given by 
(AjAj, Vi ® Vj) ioY i, j = 1, . . . ,n. The vector vi Vi is thus an eigenvector of A^A ® A^A 
of largest eigenvalue magnitude. Since vec{vk f J) = Vk®Vk we see that the supremum of 
II AX II i? is achieved for the symmetric rank one matrix vivj. From this we conclude 
that 

sup II AX A"^ Hi;' = sup IIAXA'^lli^ 

xeR"X",||x||F=i xes,\\x\\F=i 

and thus also 

sup ||(A(g) A)x|| = sup ||(Ma)x||. 

I|x||=i 11x11=1 

Notice that M^^ = (M^Jg- and so 

IKMaJ.II = IIM^JI = \\A^(S)AJ = \\{A®Ai)J 

and it then suffices to apply the definition of the joint spectral radius to conclude the 
proof. 

□ 

A matrix A is positive semi-definite (which we denote by A ^ 0) if v'^Av > for all 
f G R'^. The set of symmetric positive semi-definite matrices, denoted 8+, is a cone. It 
is not a proper cone of R,"^" because it has empty interior in R"^"; but it is a proper 
cone of the set of symmetric matrices S. If X G S+ then we clearly have AXA^ G S+ 
and so the operator X — ^ AXA^ leaves the proper cone S+ invariant. Combining this 
observation with Theorem El and Theorem |31 we deduce: 

Theorem 5. Let {Ai G R*^^"- : i = 1, . . . ,m} and denote by Ma a matrix representation 
of the linear operator X AXA^ : S ^ S. Then 

(4.3) ^ p"\Ma, + ... + MaJ < p(Ai, . . . , A„) < p'/'iMA, + ... + MaJ. 

Note that for m = 2 the bounds ()4.3|) are not so bad. Indeed, l/V^ > 0.7 and so the 
relative accuracy of the approximation is at least 70%. The main interest of this result 
resides however in the possibility of applying it recursively. Starting from an initial matrix 
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A G R"^" we define a sequence of operators acting on spaces of increasing dimensions. 
For Z = 1 we define 

and for / > 2 we define tlie operators recursively 

AV\x) = M^[,-i](x) X e Ei = S{Ei_i). 

Assume for tlie simplicity of the presentation that m = 2 and let Ai,A2 G R"^". From 
Theorem |3] we know that 

p{Af,Af)=p^\A„A,) 

On the other hand, for any / > 1 the operators Af^ and leave the cone S+(£'i_i) 
invariant. Therefore, in view of Theorem ^ we have 

|p(aw + a?)<p(4'U?')<pK + 4'0- 

Combining these last two expressions, we get the following bounds: 



l\l/2' 



(^) 



p A 



1/2' 



<p{Ai,A2) < 



p A 



1/2' 



Note that 



l-i#<(|)^/"<l 



and thus the improvement in the quality of our approximation is quite fast. Unfortunately, 
the dimensions of the spaces Ei are also growing fast; we have 

ni+i = \ni{ni + 1) 

(where rii = dimi?j) and therefore asymptotically we have ni = O (^{n/2y'^. Let us 

display, for a pair of matrices, the relative accuracy of our approximation as a function of 
the resulting dimension. 



Steps 


Accuracy 


n = 2 


n=10 


n = 100 


1 


0.707 


3 


55 


5050 


2 


0.840 


6 


1540 


* 


3 


0.917 


21 


118570 


* 


4 


0.957 


231 


* 


* 


5 


0.978 


26796 


* 


* 



We use the symbol * to mark the cases for which the dimension of the auxiliary problem 
goes beyond the abilities of modern computers. 
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5. Numerical implementation and examples 

The recursive semi-definite hfting for obtaining approximations of the joint spectral 
radius of increasing accuracy may be difficult to implement because the matrix Ma is not 
easy to express in terms of the matrix A. Consider for example the case of 2 x 2 matrices. 
Let 

ail ^12 

(3,21 (322 

then one easily compute 



A 




2aiiai2 (312 

Ma = I (311021 (3ll(322 + (3l2(321 (3l2(322 

2a21(322 (322 

In general, the matrix Ma can be expressed as follows: 

(5.1) Ma = {Q^Q)-^Q''P-\A ® A)PQ 

where P is a particular permutation matrix of size and Q is given by 



Q 



( In 


\ 





/n(n-l) 




2 

I n(n—l) , 
2 / 



(the matrix Ik is the identity matrix of size k). One difficulty with the expression (jS.lj) is 
that the permutation matrix P is tedious to construct; it essentially relates the column 
decomposition of a matrix with its diagonal decomposition. The permutation matrix 
can of course be constructed but the implementation of this construction is somewhat 
cumbersome. A way to bypass this difficulty can be achieved by using the semi-definite 
lifting only once, and then apply only Kronecker liftings. The semi-definite lifting can be 
performed either at the beginning or at the end of the process, leading to the bounds 

(5.2) p(M|; + --- + Mfj< /'(Ai, ...,A^)< p{M% + --- + Mfj 
and 

(5.3) ^ piM^^i + . . . + M^«0 < P'\^u • • • , ^m) < P(M^«. + • ■ ■ + M^g,0. 

The validity of these inequalities results from the combination of Theorem ^ Theorem 
121 Lemma 0] and Theorem EJ The expressions ()5.2|) and ()5.3|) both provide a relative 
accuracy of (l/m)-*^/*^^') but involve matrices of different size. The matrices M®' have 
size {n{n + l)/2)' whereas the matrices M.®; have size n'-{n'' + l)/2 > {n{n + l)/2)'. In 

i 

addition to this, the matrices M^' are easier to compute than M^(si since they necessitate 
the evaluation of semi-definite liftings of smaller matrices. 

Of course, other combinations of Kronecker and semi-definite liftings are also possible. 
Whenever a k-th Kronecker power is used, the current relative accuracy is k-th rooted, and 
whenever a semi-definite lifting is used, the relative accuracy is square rooted. Moreover, 
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at least one semi-definite lifting is needed if the original matrices do not leave a proper 
cone invariant. 

For the numerical implementation of these approximations notice also that the spectral 
radius of Ma^ + ■ ■ ■ + is the spectral radius of the hnear operator 

B:X-^ AiXAj + ■■■ + AmXAl^ : S ^ S 

and it can therefore be found by standard linear algebra techniques, or by solving the 
following linear inequality optimization problem 

(5.4) inf {r : rX ^ A^XAj + ■■■ + A^XA^, X y O} 

for which efficient convex optimization algorithms are available (see, e.g., ^2|) that are 
implemented in standard commercial softwares such as the Matlab "LMI Control Tool- 
box". 



6. Computational complexity analysis 

The joint spectral radius can be approximated to arbitrary accuracy. It is proved in 
that, unless P=NP, approximating algorithms of relative accuracy 1 — e for pairs of 
matrices cannot run in time polynomial in the size of the matrices and in ln(l/e). It was 
later noticed (see the note 9 on page 1260 of that a careful examination of the proof in 
[32] shows that, unless P=NP, there are no approximation algorithm of relative accuracy 
1 — e that run in time polynomial in the size of the matrices and in 1/e. This is however 
not proved in jU]. In this section we prove this result by extracting the essential part of 
the proof in [32] . 

We proceed by reduction from the classical NP-complete satisfiability problem 3-SAT. 
This problem is defined as follows. Consider a set {xi, . . .a;„} of Boolean variables. A 
literal is either a Boolean variable Xi, or its negation Xj. A three-literal clause is a disjunc- 
tion of three literals (e.g., X3 Vxs Vxe). The 3-SAT problem is the problem of determining, 
for a given collection of three-literal clauses, if there exists a truth assignment for the vari- 
ables that simultaneously satisfies all clauses. This problem is known to be NP-complete; 
finding a polynomial time algorithm for the problem is exactly equivalent to proving that 
P=NP. 

The complexity result proved in [35] is based on the construction of a pair of matrices 
which we briefiy outline (the complete construction can easily be extracted from [22] )• 
From a given instance of 3-SAT with n variables and p clauses two square matrices Ai, A2 
are constructed in polynomial time. The matrices have their entries in {0, 1}, are of size 
(n + l)(p + 1), and have a joint spectral radius that satisfies: 

p{Ai^A2) > if the instance is satisfiable 

p(Ai, A2) < (p— if the instance is not satisfiable. 
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Any algorithm of relative accuracy (1 — allows to make the distinction between 

these two cases. Moreover, since 

-^Nl/(n+2) ^ 



1 - - < 1 < 1 

V) p{n + 2) 

it is clear that an approximation algorithm of relative accuracy l — l/p{n + 2) also allows to 
decide 3-SAT. Since the size of the matrices in the construction are given by (n + l)(p + l), 
we deduce: 

Theorem 6. Unless P=NP, the problem of approximating the joint spectral radius of two 
square matrices with {0, 1} entries and with relative accuracy 1 — e cannot he obtained in 
time polynomial in the size of the matrices and in 1/e. 

From Theorem El we know that the spectral radius of two matrices with nonnegative en- 
tries is given by p(Ai, A2) = limfc^oo p^^''{Af''+Af'') and so the quantity limfc^oo p^^''{Af''+ 
Af^) is NP-hard to approximate in the sense given in Theorem IHl Theorem IHl provides 
also a rate of convergence for the approximations p^^'^^Af'' + Af''). From this we may 
prove the following complexity result. 

Theorem 7. The problem of determining, for a given integer k > and for a given pair 
of matrices Ai,A2 with nonnegative entries, if 

p{Af + Af) < 1 

is a problem that is NP-hard. 

Proof. We proceed by reduction from 3-SAT. From a given instance of 3-SAT with n 
variables and p clauses we use the construction given in [32] to construct two square 
matrices Bi,B2 that have their entries in {0, 1}, are of size (n -|- l){p + 1), and have a 
joint spectral radius that satisfies: 

p{Bi,B2) > pV(n+2) iftheinst ance is satisfiable 

p(Bi, B2) < (p — if the instance is not satisfiable. 

We then choose a G R with 0<a<l, rGQ and k >0 such that a < (1/2)^/'^' and 

1 1 1 1 1 

(6.1) (p — l)"+2 < — (p — l)"+2 < r < Q;g"+2 < p"+2. 

a 

Consider now the matrices Ai = (l/r)Bi and A2 = (l/r)i?2. We claim that p(Af^ -|- 
Af'^) < 1 iff the instance of 3-SAT is not satisfiable. Indeed, assume first that the 
instance of 3-SAT is not satisfiable. Then 

p{B„B2)<{p-lY/^''^'\ 

Since p{Af^ + Af^) is an approximation of relative accuracy 1/2^^ of p{Ai, A2), we have 

p{Af + Af) < 2^V(^i, A2) < —p{Bi, B2) <—ip- 
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and by using the inequalities (jfi.lj) we conclude p{Af'^ + ^) < 1. In an analogous way 
one see that, if the instance of 3-SAT is satisfiable, then p^Af'' + > 1 and the proof 
is therefore complete. 

□ 



7. Ellipsoid approximation 

Let us now compare our approximation of the joint spectral radius with another ap- 
proximation appearing in the literature. The following approximation, called ellipsoid 
approximation, is introduced in P: 

.1/2 



(7.1) p{A,,...,Ar, 



inf {r : tX >z AiXAj , i = 1, . . . , m, X ^ 0} 



X 



This expression corresponds to the best ellipsoid norm for the set of matrices. In jB] the 
following inequalities are proved 

^ p(Ai, ...,Am)< p(Ai, . . . , A^) < p(Ai, ...,Am). 

Thus, the ellipsoid approximation has a relative accuracy of 1/ i/n. This accuracy de- 
creases with the size of the matrices but does not depend on the number of matrices 
in the set. Using our results we can prove that the relative accuracy of the ellipsoid 
approximation is in fact bounded by 1/ -/m where m is the number of matrices. 

Theorem 8. Let {Ai : z = 1, . . . , m} he a set of matrices and define the ellipsoid approx- 
imation p by \7.1\ ). Then 

^ p{Ai, ...,Am)< p{A,, ...,A^)< p{Ai, ...,Am). 

Proof. It is easy to see that p{Ai, . . . , Am) < p(^i, • • • , Am). For proving the first inequal- 
ity, note that the spectral radius of the linear operator 

m 

can be represented as follows: 

m 

(7.2) p{B) = inf {r : rX ^ J] AXAJ, X y 0}. 

i=l 

If a pair (X, r) is feasible for the optimization problem in (j7.2j) . then it is also feasible for 
the optimization problem in (j7.ip . Therefore 

p{A„...,AJ<p'/\B). 

Finally, by Theorem El we have 

^pi/2(S)<p(Ai,...,/l^) 
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and so we get the following bounds 

(7.3) ^/5(Ai, . . . , A^) < ^p"\B) < p(Ai, . . . , A™) < p{A,, ...,A^)< p^l\B). 

Thus we see that the ellipsoid estimate . . . , Am) has relative quality 1/ y/m. □ 
By using Theorem |21 one can obtain a statement analogous to that of Theorem El 

Corollary 3. Let {Ai : i = 1, . . . , m} he a set of matrices and define the ellipsoid approx- 
imation p by \7.1\ ). Then 

(7.4) ^ p"\Al ...,Al)< p(Ai, . . . , A^) < p''\Al ...,Al). 

For — 1 , this result is also proved in jTj . For arbitrary k a proof can be reconstructed 
from the result proved in [21]. Notice however that approximation algorithms that are 
directly based on (|7.4j) are not efficient because they require to solve linear matrix in- 
equalities of very large dimension. Approximations based on (j5.2j) are easier to obtain 
because there exists many efficient numerical algorithms for approximating the spectral 
radius of a matrix. 

8. Discussion 

The results presented in this paper introduce the possibility to compute approximations 
for the joint spectral radius with worst-case theoretical guarantees. Of course, the com- 
putational complexity of these estimates needs further examination. Indeed, the spectral 
structure of the matrix A®^ is very simple and this simplicity must be inherited somehow 
by the matrix Af^ + Af^. It therefore appears to be an interesting problem to investigate 
the possibility of constructing efficient procedures for finding the spectral radius of this 
sum; of course there are theoretical limitations on what can be achieved since we have 
shown that the problem of deciding p{Af^ + Af") < 1 is NP-hard. 

Another interesting issue is that of determining if worst-case theoretical guarantees can 
be provided for the approximation 

max p{A„f'^<p{Ai,...,Am). 

a&{l,...,m}'^ 

In many numerical examples, this approximation does in fact perform at least as well. 
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